Testing rig for new parallel solver

In [1]:

    

from zephyr.Frontend import *


    

In [2]:

    

import matplotlib.pyplot as plt
import matplotlib.cm as cm
%matplotlib inline

import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 150 # Change this to adjust figure size

# Plotting options
font = {
    'family': 'Bitstream Vera Sans',
    'weight': 'normal',
    'size': 8,
}

matplotlib.rc('font', **font)


    

In [3]:

    

cellSize    = 1             # m
freqs       = [2e2]         # Hz
density     = 2700          # units of density
Q           = np.inf        # can be inf
nx          = 164           # count
nz          = 264           # count
freeSurf    = [False, False, False, False] # t r b l
dims        = (nx,nz)       # tuple
nPML        = 32
rho         = np.fliplr(np.ones(dims) * density)
nfreq       = len(freqs)    # number of frequencies
nky         = 48            # number of y-directional plane-wave components
nsp         = nfreq * nky   # total number of 2D subproblems

velocity    = 2500          # m/s
vanom       = 500           # m/s
cPert       = np.zeros(dims)
cPert[(nx/2)-20:(nx/2)+20,(nz/2)-20:(nz/2)+20] = vanom
c           = np.fliplr(np.ones(dims) * velocity)
cFlat       = c
c          += np.fliplr(cPert)
cTrue       = c

srcs        = np.array([np.ones(101)*32, np.zeros(101), np.linspace(32, 232, 101)]).T
recs        = np.array([np.ones(101)*132, np.zeros(101), np.linspace(32, 232, 101)]).T
nsrc        = len(srcs)
nrec        = len(recs)
recmode     = 'fixed'

geom        = {
    'src':  srcs,
    'rec':  recs,
    'mode': 'fixed',
}

cache       = True
cacheDir    = '.'

# Base configuration for all subproblems
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'c':        c.T,        # m/s
    'rho':      rho.T,      # density
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
    'geom':     geom,
    'cache':    cache,
    'cacheDir': cacheDir,
}

# Generate the configuration update objects, which contain freq, ky and weight
subConfigSettings = {
    'freqs':    freqs,
    'nky':      nky,
    'cmin':     c.min(),
}

handles25D = getHandles(systemConfig, subConfigSettings)
forward25D = handles25D['forward']
clear25D = handles25D['clear']


    

In [4]:

    

%%time
fields = []
data = []
for isrc in xrange(101):
    lfield, ldat = forward25D(isrc, False)
    fields.append(lfield)
    data.append(ldat)
data = np.array(data)

dTrue = data
fTrue = fields


    

    




CPU times: user 39.7 s, sys: 11.3 s, total: 51.1 s
Wall time: 6min 23s

In [5]:

    

clip = max(abs(fields[0])) * 1e-1

fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
plt.imshow(fTrue[0].reshape((nz+1,nx+1)).real, cmap=cm.bwr, vmin=-clip, vmax=clip)
ax1.invert_yaxis()
ax2 = fig.add_subplot(1,2,2)
plt.imshow(fTrue[-1].reshape((nz+1,nx+1)).real, cmap=cm.bwr, vmin=-clip, vmax=clip)
ax2.invert_yaxis()


    

In [6]:

    

fig = plt.figure()
ax1 = fig.add_subplot(2,2,1, aspect=1)
plt.imshow(dTrue.real, cmap=cm.bwr)

ax2 = fig.add_subplot(2,2,2, aspect=1)
plt.imshow(dTrue.imag, cmap=cm.bwr)

ax3 = fig.add_subplot(2,2,3, aspect=1)
plt.imshow(abs(dTrue), cmap=cm.jet)

ax4 = fig.add_subplot(2,2,4, aspect=1)
plt.imshow(np.angle(dTrue), cmap=cm.gray)

fig.tight_layout()


    

In [7]:

    

systemConfig.update({'c': cFlat})

handles25D = getHandles(systemConfig, subConfigSettings)
forward25D = handles25D['forward']
gradient25D = handles25D['gradient']
clear25D = handles25D['clear']


    

In [8]:

    

%%time
#fields = []
data = []
for isrc in xrange(nsrc):
    #lfield, ldat = forward25D(isrc)
    ldat = forward25D(isrc)
    #fields.append(lfield)
    data.append(ldat)
data = np.array(data)

# source receiver
dFlat = data
#fFlat = fields
dResid = dFlat - dTrue


    

    




CPU times: user 36.3 s, sys: 7.1 s, total: 43.4 s
Wall time: 6min 11s

In [9]:

    

%%time
grad = []
for isrc in xrange(nsrc):
    dR = dResid[isrc,:]
    grad.append(gradient25D(isrc, dR))
grad = np.array(grad).reshape((nsrc, nz+1, nx+1))
gradSum = grad.sum(axis=0)


    

    




CPU times: user 1min 10s, sys: 18.6 s, total: 1min 29s
Wall time: 11min 37s

In [10]:

    

clipper = abs(gradSum).max()

fig = plt.figure()

ax1 = fig.add_subplot(1,3,1, aspect=1)
plt.imshow(c.T, cmap=cm.jet, vmin=velocity, vmax=velocity+vanom)
plt.title('True model')

ax2 = fig.add_subplot(1,3,2, aspect=1)
plt.imshow(gradSum.real, cmap=cm.bwr, vmin=-clipper, vmax=clipper)
plt.title('Full gradient')

ax3 = fig.add_subplot(1,3,3, aspect=1)
plt.imshow(grad[0].real, cmap=cm.bwr, vmin=-clipper/nsrc, vmax=clipper/nsrc)
plt.title('Single-source gradient')

fig.tight_layout()


    

In [11]:

    

clear25D()


    

    




Cleared stored matrix terms for 48 systems.

In [11]: